Thermalisation of Local Observables in Small Hubbard Lattices 
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We present a study of thermalisation of a small isolated Hubbard lattice cluster prepared in a pure 
state with a well-defined energy. We examine how a two-site subsystem of the lattice thermalises with 
the rest of the system as its environment. We explore numerically the existence of thermalisation 
over a range of system parameters, such as the interaction strength, system size and the strength of 
the coupling between the subsystem and the rest of the lattice. We find thermalisation over a wide 
range of parameters and that interactions are crucial for efficient thermalisation of small systems. 
We relate this thermalisation behaviour to the eigenstate thermalisation hypothesis and quantify 
numerically the extent to which eigenstate thermalisation holds. We also verify our numerical results 
theoretically with the help of previously established results from random matrix theory for the local 
density of states, particularly the finite-size scaling for the onset of thermalisation. 
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I. INTRODUCTION 

Understanding the quantum origins of statistical me- 
chanics has seen renewed interest over the last few years, 
in part motivated by experimental progress in degenerate 
atomic gases, but also due to independent theoretical ad- 
vances The central question is as follows. Consider 
a closed quantum system prepared in a pure quantum 
state. Does it evolve in time to a thermal state? If so, in 
what sense is it a thermal state? 

In this paper, we focus on observables that are local 
to a subsystem of the full system. Thus, we discuss the 
'thermalisation' of this subsystem with the rest of the sys- 
tem as a bath (see Fig. [T|) . We will discuss the conditions 
for the eventual thermalisation of this subsystem. This 
has been studied in many systems and we will study 
a system of interacting fermions in this context. The pic- 
ture of a local subsystem in a closed system also naturally 
makes contact with the conventional framework of statis- 
tical mechanics where thermal equilibrium is achieved by 
a weak coupling AV^ between a system and its environ- 
ment. 
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FIG. 1. Schematic diagram of a closed system divided con- 
ceptually into a subsystem and a bath. 

Thermalisation in closed quantum systems has been 
shown to have its origins in entanglement. To be specific, 
let us consider a composite system with a Hamiltonian 
H = Hs + Hb + W where Hs,b describes the dynam- 
ics of subsystem {S) and the bath {B) respectively while 
AV^ couples the subsystem to the bath. The exact eigen- 



states \A) of this Hamiltonian are typically superposi- 
tions of many eigenstates of the decoupled system (A = 0) 
which are products of the subsystem and bath states: 
1^) = y^jf. Csb|s)<; ® \b)B- The idea of 'canonical typical- 
ity' |lOl4l4l | states that almost any pure state composed 
of many energy eigenstates \A) within a narrow energy 
window will give rise to a canonical distribution for the 
measurements of local or few-body observables within the 
subsystem. This emerges because the pure state is an en- 
tangled combination of subsystem and bath eigenstates. 
We will consider in this paper a system prepared initially 
at time i = in a pure state that is a product state of 
the subsystem and bath states. Such a state is typi- 
cally a superposition of many closed-system eigenstates: 
\^{t = Q)) = \(l,)s(E)\i^)B=EAdA{t = 0)\A). While this 
initial state is special and cannot be considered as 'typ- 
ical', we expect the wavefunction will, in general, evolve 



in time (|5'(i)) = 



^(0))) towards a state that falls 



into the domain where canonical typicality applies. The 
sufficient conditions for this to occur have discussed in 
recent papers In this paper, we will investigate 

conditions for thermalisation in a small Hubbard-model 
system. 

An alternative view is the eigenstate thermalisation 
hypothesis [iMl (ETH). The time evolution of any 
few-body observable {^{t)\0\'^{t)) involves the inter- 
ference of eigenstates at different frequencies: (O) = 
J2ab d*AdB{A\0\B)e''^^^-^^'>* where Ea,b is the energy 
of the eigenstates \A) and \B). The eigenstate thermal- 
isation hypothesis says that destructive interference re- 
moves all A ^ B terms and that 



{A\0\A) « (O) 



Ea 



(1) 



where the right-hand side denotes the thermal average of 
O when the total system has energy Ea- This paints a 
very different picture of thermalisation compared to the 
scenario for classical statistical mechanics where states 
diffuse ergodically through phase space constrained by 
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energy conservation. 

These concepts are powerful because they guarantee 
thermahsation for closed quantum systems. They de- 
pend crucially on the very high dimensionality of the 
Hilbert space of quantum states. In this paper, we aim to 
gain insight into these ideas by testing the limits of these 
hypotheses in terms of the breakdown of thermahsation 
for a small closed quantum system. We take our moti- 
vation from cold atom experiments with optical lattices 
and single-site addressibility [13, HH . We choose a lattice 
system of interacting fermions in a normal metallic state. 
In particular, we present a study of the thermahsation of 
a composite system consisting of a 2-site subsystem and 
a {L ~ 2)-site bath in a one-dimensional Hubbard ring. 
We avoid the issue of integrability and the generalised 
Gibbs ensemble [l^ by choosing parameters such that it 
is a non-integrable system. 

In order to study the thermahsation of the subsystem, 
we need to calculate the long-time behaviour of reduced 
density matrix p of the subsystem: 

p{t)^TrB\^{t))(^it)\. (2) 

where Tr b denotes a trace over the bath degrees of free- 
dom. A thermalised system corresponds to a diagonal 
reduced density matrix with diagonal elements given by 
the Gibbs distribution. We explore whether thermahsa- 
tion occurs over a range of system parameters. We find 
numerically (section IIIip that thermahsation occurs in 
surprisingly small systems. For a system of a given size, 
there is a threshold for the onset of thermahsation, both 
in terms of the coupling strength A and the interaction 
strength. In particular, we study the size dependence 
of the threshold Ath that the coupling strength has to 
exceed to achieve thermahsation. We demonstrate that 
this threshold for thermahsation agrees with the ETH 
criterion Q for thermahsation (sections IIII G[) . Indeed, a 
theoretical threshold Aeth determined from the ETH cri- 
terion has the same size dependence as the empirical Ath 
(section lV Ap . We also argue that both of these thresholds 
mark the onset of non-perturbative mixing of eigenstates 
due to the subsystem-bath coupling (at a threshold Anp). 

From a separate perspective we can study 

the thermahsation process in terms of the statistics of 
the eigenstates. We study the statistics of the over- 
lap {A\sb) of the eigenstates \A) of the coupled system 
with the eigenstates of the decoupled system which are 
product states |s6) = |s)s €5 Interestingly, at weak 

subsystem-bath coupling where the onset of thermahsa- 
tion occurs, the distribution for the overlaps fits a hyper- 
bolic secant distribution (section IIVAI). T his is in con- 
trast to previous conjectures fvA llSl |23| from random 
matrix theory which suggest that these types of overlaps 
should follow a Normal distribution at weak coupling. 

Using our results for the overlap distribution, we show 
numerically (section IIII Gp that the eigenstate thermal- 
isation hypothesis ([l} holds for the projection operator 
Ps = J2b which projects onto the subsystem state 

s in the parameter regime where the subsystem is ther- 



malised. This can be explained theoretically (section 
lIVBp using known results for the variance of the overlap 
distribution. This observation for P, then leads directly 
to a thermalised reduced density matrix. (See section 

mm) 

In this paper, we also highlight the importance to ther- 
mahsation of the strength of interaction within the bath. 
We find that, at least for our small bath and subsystem, 
a finite interaction strength is needed for thermahsation. 
This is consistent with the expectation that thermahsa- 
tion is aided by inelastic scattering in the bath. 

We point out that, although we have focussed on the 
thermahsation of a spatially local subsystem, one can also 
study the thermahsation of few-body observables over 
the entire system. This has been studied particularly 
in the context of quantum quenches in a variety of sys- 
tems p^[24M27j . including integrable systems (23 . [28l - l33| . 
Moreover, one can discuss the d yna mics of the relaxation 
towards a thermal state [13, [Ull^] • Both of these issues 
are beyond the scope of this paper. 

This paper proceeds as follows. The following section 
introduces the Hubbard model we study, discusses how 
the system is prepared initially and provides a framework 
for studying thermahsation. In Section IIIIl we present 
a comprehensive set of results for the thermahsation of 
two-site subsystems in small Hubbard rings. We con- 
sider the effects of subsystem-bath coupling strength on 
thermahsation and link these results to eigenstate ther- 
mahsation. We further demonstrate the role of inter- 
actions between fermions before exploring system size 
dependence and, finally, the energy width of the initial 
prepared state. Section IIVI introduces results from ran- 
dom matrix theory concerning the nature of the eigen- 
states of the coupled system. From these, we review the 
arguments leading to eigenstate thermahsation [l^ and 
derive the scaling behaviour associated with the close- 
ness to perfect eigenstate thermahsation. In Section fVl 
we present an account of system-size scaling by consider- 
ing a threshold for non-perturbative mixing of uncoupled 
composite eigenstates for the cases of both interacting 
and non-interacting fermions. In Section IVIl we discuss 
implications for experiments. Finally, in Section [VIIl we 
give our conclusions. 



II. THE MODEL AND THEORETICAL 
FRAMEWORK 

We will consider closed quantum systems with unitary 
time evolution. The system is prepared in an initial state 
|5'(t = 0)) which evolves in time |^'(t)) = e"*-^*|*(0)) 
under the influence of the Hamiltonian H. In this sec- 
tion, we discuss the specific model studied in this work 
and our choice of initial states. Since we will investigate 
thermahsation of a subsystem of this closed system, we 
will also discuss our criteria for a thermal state. 
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A. Hubbard Hamiltonian 

We divide the system into a local subsystem (5) and a 
bath (B). The subsystem (bath) is described by a Hamil- 
tonian Hg (Hb) acting on the subsystem (bath) Hilbert 
space. Let us denote the subsystem (bath) eigenstates as 
\s)s ( 1 6) b) with energies Eg (cb). The subsystem and bath 
are coupled by a Hamiltonian AT^ . We will use A as a tun- 
able parameter to control the strength of this coupling. 
At A = 0, the eigenstates are products of subsystem and 
bath eigenstates, \sb), with energies Egb = + At 
non-zero A, the eigenstates are in general entangled with 
respect to the subsystem-bath partition. We denote these 
composite eigenstates by \A) (using an uppercase index) 
and their energies by Ea- 




Subsystem 



FIG. 2. Schematic diagram of a two-site subsystem in a lattice 
with 9 sites (i = 9). 

In this work, we focus on the Hubbard model away 
from half filling as a simple model of interacting fermions. 
More specifically, we consist of a two-site subsystem in 
an L-site Hubbard ring of fermions such that the Hamil- 
tonian takes the form H = Hg + Hb + XV with 

-ffs = - ^ Ja{cic,C2cr + h.c.) + U{nifnii + n2fn2i) , 
Hb = - ^ Jc{c\„Ci+i.a + h.c. 

1 = 3 (T=t,i 
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FIG. 3. The density of states g{Eo) of the system at composite 
energy Eo for different coupling strengths, A (as labelled), for 
an L = 9 site lattice where U = J = 1. g{Eo) is generated 
as a histogram by counting eigenstates in a Gaussian window 
centered on Eq with width 0.5J. 
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FIG. 4. The density of composite states g{Eo) of the sys- 
tem at composite energy Eo for different Hubbard interaction 
strengths, U (as labelled), for an L = 9 site lattice where 
A — 0.5. g{Eo) is generated as a histogram by counting eigen- 
states in a Gaussian window centred on Eq with width 0.5. 
(J=l.) 



where c-^ is a creation operator for a fermion with spin 
a at site i and riia- — cl^Cia is the number operator on 
site i with spin a. This Hamiltonian describes a ring 
with the subsystem sites i = 1,2 and bath sites i = 3 to 
L with two links between the subsystem and the bath. 
Note that, in the case of A = 1, the Hamiltonian describes 
an homogeneous ring. We choose the hopping integrals 
Ja = J{1 + £, sgn(cr)), with ^ — 0.05 to remove level de- 
generacies associated with spin rotation symmetry. (We 
will use J as the unit of energy.) Breaking spin sym- 
metry and the presence, in general, of modified hopping 
integrals between sites i = 2 and 3, as well as between 
sites i = L and I make this system non-integrable for 
non-zero U. 



The total particle number, N, and spin component, 
S"^ , are conserved in addition to the total energy of the 
composite system. In the numerical results we present, 
we consider lattices with up to L = 9 sites and with 8 
fermions of total spin = 0. The two-site subsystem 
has A/5 — 16 eigenstates and the 7-site bath has 8281 
eigenstates, while the composite 9-site system has a total 
of AI — 15876 states and an average level spacing A ~ 
10-3 J. 

The spectrum of the composite system has a smooth 
quasi-continuous density of states giE^) for a range of A 
and U. This is illustrated in Figs. [3] and 21 The centre 
of the spectrum is located at Eq ~ 1.77 J. The spectrum 
develops peaks for large A and U/ J. In the case of A ^ 1, 
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we attribute this to single-particle states with a large en- 
ergy splitting proportional A J on the two links connect- 
ing the subsystem and bath. For [//J ^ 1, we attribute 
the peaks to a large energy gap to doubly-occupied sites 
(sometimes referred to an 'upper Hubbard band' in the 
theory of strongly electron systems or doublons in the 
cold atoms literature). 



B. Initial States 

Throughout this work, we consider the composite sys- 
tem to be prepared in a pure state which is a product 
state of a subsystem state and a bath state: 

\^{t^Q),Eo)^\ct>)s®Y. (4) 

The initial subsystem state, \4>)s, can be, for example, 
|t,i)s which is prepared with antiparallel spins on sites 
1 = 1 and 2 of the lattice. The initial bath state contains 
a linear combination of B bath eigenstates \hi)B- These 
bath states are chosen to be within an energy window of 
bath states such that = Eq. The energy window 

has a width in energy of (5b (specified by the state indices 
in the range 6/ < 6 < &„). Unless stated otherwise, we 
choose 5b — 0.5 J which is small on the scale of variations 
in the density of states. For a 7-site bath, this window 
contains about 100 bath eigenstates. 

As the system evolves in time, the state of the sub- 
system can be described by the reduced density matrix 
(RDM) as defined in equation We will study the 
reduced density matrix elements using the subsystem 
eigenstates \s)s as the basis: p^s' = s(s|p|s')s- We ob- 
tain the wavefunction |^(i)) of the composite system us- 
ing the eigenstates and energy eigenvalues from the exact 
diagonalisation of the Hamiltonian H: 

|vI'(t))=^e-^-*|A)(A|vI/(0)). (5) 

A 



C. Thermalisation 

To assess whether the subsystem reaches a thermal 
state, we must be more precise about the criteria for a 
thermal state. We start with the conventional definition 
of thermal equilibrium in the canonical ensemble. To de- 
fine this 'canonical thermal state', oj, of the subsystem, 
we set up the composite system at a total energy Eq in 
a microcanonical mixed state and consider the regime 
where the coupling between the subsystem and the bath 
is negligible. In this way, the thermal state of the sub- 
system may be determined by counting bath states in 
an energy window, conserving the total energy and the 
global and particle number N . The reduced density 



matrix is diagonal and is given by Wss = s{s\'-^\s)s' 

^ Mb{Eo^es,N-n,,S^--sl) 

Es'Mb{Eo-es',N-ns>,S^-sl,) ^' 

where Mb{eb, nbjSl) is the number of bath states with Ub 
fermions of spin in a window of width ~ Sb centred 
on energy e^. This thermal RDM is a function of e^, Hs 
and Sg , arising from the global conservation laws of the 
system. For fixed and and for a range of subsystem 
energies small compared with features in the density 
of states, the smooth density of states allows us to write 
the RDM in the Boltzmann form uigs ~ e^'^"^'^ , with the 
inverse temperature given by 

1 ^ dlogMBieb,N~n,,S^ ^st) 
T deb 

We can in principle deduce a chemical potential and Zee- 
man field by considering variations in Us and s^. How- 
ever, we are considering small systems where the discrete- 
ness of these quantities cannot be ignored and Mb is not 
a smooth distribution of and Ug. Nevertheless, Eq. ([6]) 
may be used to specify a thermal state of the subsystem 
for, in principle, any bath size. 

We note that such a Gibbs-like distribution has just 
three parameters, differing from the 'generalised Gibbs 
distribution' for integrable systems [l^, HI, HI], where 
the number of parameters extends with system size. 

Let us now turn to the RDM that we obtain from the 
unitary evolution from an initial pure state. We will be 
examining the behaviour of the RDM at long times. It 
is useful to define the time average: 

rss - lim - / dt' p{t') (8) 

t >00 t Jq 

If the reduced density matrix reaches a steady state at 
long times, this state will be equal to the time average 
r. We expect this to become diagonal. Using Eq. ©, we 
see that the diagonal elements of the RDM are given by 

{s\p\s) ^ J2 e-''^^'''^^^Hsb\A){A\^{0)){'^{0)\B){B\sb) . 

ABb 

(9) 

Averaging over time for long times identifies Ea with 
Eb- Since we have lifted all symmetry-related degenera- 
cies, this also identifies states A and B in the sum above 
(barring accidental degeneracies). So we see that 

r,,=^|(*(0)|^)HA|P,|^) 

A 

with = XI ■ <^^0) 

b 

The operator Pg projects from the composite Hilbert 
space on to the subsystem state \s)s by tracing over bath 
states. 

The coefficients (^'(0)|yl) contain the information 
about the initial state. However, this steady state may 
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still be very close to a state which is independent of 
initial conditions. A sufficient ondition is given by the 
'eigenstate thermalisation hypothesis' [l^. Recall that 
from Eq. Q our initial state |5'(0)) has been set up 
within a narrow energy window. So, the overlap of 
(^'(0) I A) should be only non-zero in a window of eigenen- 
ergies. (We will give a more quantitative discussion of 
the width of this window in Section IIV Al ) The eigen- 
state thermalisation hypothesis assumes that (AjPslyl) 
for a system with composite energy Eq depends only 
weakly on the choice of \A) in this window of eigenen- 
ergies. This allows us to replace (A|Ps|^) by its average 
value (^|Ps|^) over the eigenenergy window. In that 
case, r,, ~ (A|P,|A) I (*(0)l^) I' = {A\PM)- Thus, 
we see that the steady state r^s may indeed be indepen- 
dent of initial conditions. 



Furthermore, if (A|Ps|A) is close to the value uJssiEo) 
for the canonical ensemble ©, then Eq. ([TU)) can be writ- 
ten as — ^ssiEo). Thus, will be close to the canon- 
ical thermal state to for any initial state with a definite 
energy. 

In summary, we break down the question of whether 
the subsystem thermalises into four criteria similar to the 
ones in Ref. [iJI. Our criteria are 



1. Firstly, we should establish that the reduced den- 
sity matrix reaches a steady state at long times. 



2. The steady state should be diagonal in the subsys- 
tem energy eigenbasis with all off-diagonal elements 
falling to zero for long times. This demonstrates a 
loss of quantum coherence. 
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FIG. 5. Left: Time dependence of the initial state occupa- 
tion probability pss for initial state \s)s = | ti t) for three 
coupling strengths A = 0.05, 0.5 and 0.3. Right: a measure 
of the magnitude of the off-diagonal elements Eod, defined 
by Eq. Interaction strength U — J, width of the initial 

bath state 5b = 0.5J, total system energy Eo — — 2J, size 
L = 9. 



III. NUMERICAL RESULTS 



A. Long time behaviour 



3. This steady state should have no memory of the 
initial state, such as the precise way in which the 
subsystem or the bath is prepared. 

4. Finally, we ask if this steady state is close to the 
canonical thermal state lu. If this is the case, we will 
say the system exhibits 'canonical thermalisation'. 

We are leaving open the possibility that the subsys- 
tem reaches a steady state with no memory of the initial 
state, but does not resemble the canonical thermal state. 
This may be possible since the canonical state has been 
derived assuming the bath states are unperturbed by the 
coupling with the subsystem which may not hold in the 
small quantum systems studied here when the coupling 
A is of order unity. 

We explore these questions with numerical studies of 
the Hubbard Hamiltonian in the following section. We 
will discuss the more sophisticated picture of the eigen- 
state thermalisation hypothesis separately in Section HVl 



We proceed to demonstrate that the first two criteria 
for thermalisation listed at the end of Section III Cp are 
met for a range of system parameters. These require- 
ments are that the reduced density matrix p should ap- 
proach a steady state at long times, with its off-diagonal 
elements falling to zero. 

Initial states of the form in Eq. Q were constructed 
with the initial subsystem state | tj t)s- It was found that 
evolving p{t) in time under the Hamiltonian H results in 
almost steady states for a wide range of A, provided the 
composite energy is not close to the edge of the spec- 
trum. With interaction strength U = J and bath width 
6b = 0.5 J, this range is 0.05 < A < 3 for a composite en- 
ergy Eq — —2 J. Examples are shown in the left panel in 
Fig. [5l This shows the relaxation of the diagonal element 
of p corresponding to the initial-state occupation proba- 
bility for couplings A — 0.05, 0.5 and 3. The dynamics 
of the relaxation is fast and featureless for A up to 1. 
The temporal fluctuations around the long-time steady 
state are small at this energy Eq — —2J. In fact, if we 
use an energy close to the centre of the spectrum of the 
composite system (Eq = 1.77 J), temporal fluctuations 
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are even smaller for a given A. On the other hand, for an 
energy closer to the edges of the spectrum, the density of 
states is small so that few composite states construct the 
initial state. The presence of only a few frequencies in 
the time evolution limits the closeness to a steady state 
achievable. 

We see beating oscillations at A = 3 which we at- 
tribute to the strongly split single-particle states at the 
two subsystem-bath links at large A. Indeed, for even 
larger A, p{t) no longer reaches a steady state, and the 
frequency spectrum begins to show peaks at frequencies 
which are integer multiples of A J. 

Let us now investigate the second condition for ther- 
malisation that off-diagonal elements fall to zero with 
only small temporal fluctuations, as predicted in [l^. We 
compute the root-mean-square sum of these off-diagonal 
elements: 



SodW= /El/'-'WP- (11) 

y s<s' 

This is shown in the right panel of Fig. [5] By construc- 
tion, T,oD is larger than any single off-diagonal element. 
We see that the effect of off-diagonal elements can be ne- 
glected at long times: "Sod is, with decreasing A, shown 
to be from ^ 10^^ down to ^ 10^^ times smaller than 
each diagonal element. 

Having established that, within a range of coupling 
strengths, the subsystem RDM does reach a diagonal 
steady state with only small temporal fluctuations, we 
will now use Eq. pUj) to compute the steady-state form r 
without explicitly computing p(t) at many times and tak- 
ing a time average. This is less computationally expen- 
sive and provides a definitive long-time subsystem state 
without the need for numerically averaging out small 
temporal fluctuations. 

We will now proceed to explore two further require- 
ments of thermalisation: these are the extent of initial- 
state independence and closeness to the thermal state oj. 
The effective temperature of the subsystem may also be 
estimated from r. 



B. Quantifying Thermalisation 

Next we develop measures to characterise the extent to 
which the third and fourth of our thermalisation criteria, 
listed in Section Hi CI are met. 

Criterion 3 in Section III CI is concerned with the loss 
of memory of the initial state at long times. To quantify 
the variation in the steady state due to different initial 
states, we introduce the measure Ar which measures the 
root-mean-square variation in diagonal reduced density 
matrix elements for different initial subsystem states: 

^r = lY.^{rss')-{rssff^ (12) 



with (. . .) denoting an average over all 16 initial states 
in the subsystem Fock basis, as in Eq. (IT^ . We expect 
Ar to be small when the long-time steady state no longer 
depends on how the system was initially prepared. 

Criterion 4 in Section lll Cl addresses the closeness of the 
subsystem state at long times to the canonical thermal 
state uj. We quantify this with the quantity (7^^, defined 
as: 

<^u^ = \ ^{\rss-i^ss\) ■ (13) 

s 

Here, (...) denotes an average over all 16 initial states 
in the subsystem Fock basis (eigenstates at J = 0). As 
such, this is a measure of the average distance to the 
thermal state, u, for the set of initial subsystem states 
spins localised on the lattice sites. It is a spec ial case o f 
a more general distance measure [3], (^Tr-y/ (r — w)^), 
which equals a^j in the case where the elements of r in 
the subsystem eigenbasis form diagonal matrices. As es- 
tablished above, this is the case for 0.05 < A < 3. Within 
this range, we may interpret cr^ as the probability, upon 
making measurements on the subsystem, that rss could 
be distinguished from ojss [l^ . 

From the definitions of these two measures, it is clear 
that if Ar is large then (T(j is necessarily large too: if there 
is a large variation in r for different initial states, many 
of these states must be far from the uniquely defined 
canonical thermal state uj. Conversely, it is possible for 
aui to be large with Ar small, because the subsystem may 
relax consistently to a state r other than the canonical 
state UJ. 

We will also compute the von Neumann entropy of the 
subsystem. Because off-diagonal elements of p{t) are vir- 
tually zero at long times even for very small A, we intro- 
duce an initial-state-averaged subsystem entropy for the 
equilibrium state, which we define by 

S=-J2{rss\ogrss} (14) 

s 

where (...) denotes an average over all initial subsystem 
states in the subsystem Fock basis. 

We would also like to characterise subsystems showing 
thermalisation with an effective temperature. As dis- 
cussed in Section III CI if we consider the subsystem at 
a given particle number Us and spin s^, we expect the 
steady-state RDM, r, to approach the Boltzmann form 
([7|) for u! if the subsystem relaxes to the canonical thermal 
state UJ. Therefore, we extract an effective temperature 
Tcff from the RDM, r, of the steady states that we find 
using a least-squares fit to the form 

log r ss = + const . (15) 

We will focus on the four-state subsector with rig = 2 and 
= because it is the subsector with the largest number 
of bath states. Note that it is possible that we can have a 
good fit to this form with an effective temperature even 
if the steady state is not close to the canonical state uj. 
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C. The Role of Coupling Strength 




0.01 



(S) of the subsystem. We will now discuss how these 
measures of thermalisation change over a broad range of 
subsystem-bath coupling strengths A. We show results 
at different total energies Eo between — 4J and 1.77J. 

In Fig. [51 we present our results for Ar, and S as 
a function of the coupling strength A (for a system with 
U = J and an initial state of bath width Sb = 0.5 J). Our 
results for Ar demonstrate that the subsystem reaches a 
steady state with little dependence on initial conditions 
over a wide range of coupling strengths A. We see sig- 
nificant dependence on initial state beyond this range, 
at both small and large A. Moreover, our results for cr^ 
show that the subsystem reaches the canonical state uj 
over a similar, albeit slightly narrower, range of coupling 
strengths. Outside this range, the long-time steady state 
shows strong deviation from the canonical state. 



Canonical thermal 




0.01 




0.01 



FIG. 6. Memory of initial state Ar, closeness to the ther- 
mal state a^i and the subsystem entropy S as a function 
of coupling strength A for different composite energies Eo- 
({/ = J = 1, <5s = 0.5, L = 9.) 



In the previous section, we discussed how we measure 
the memory of the initial conditions (Ar) in the steady 
state, closeness (au) to the canonical thermal state, the 
effective subsystem temperature (Teg) and the entropy 



— I 1— 
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FIG. 7. A schematic diagram indicating the range in A where 
the subsystem reduced density matrix, r, is diagonal, and 
where it is close to the canonical thermal reduced density 
matrix, uj. 

In the coupling range where Ar and are both small, 
we find that the entropy S reaches a plateau as a func- 
tion of A. Beyond this range at low A, the subsystem 
entropy S drops with decreasing A. This is consistent 
with the subsystem retaining information of its initial 
conditions. On the other hand, the entropy rises when 
A is increased beyond the plateau. The asymmetry be- 
tween low and high coupling indicates that the departure 
from thermalisation at small and large A have different 
physical origins, as we discuss later. The behaviour of 
the subsystem reduced density matrix as a function of A 
is summarised by the schematic diagram in Fig. [7] for the 
range of energies shown in Fig. [6] 

We show in Fig. |8] the effective temperature Tcff ex- 
tracted at different energies Eq using the fit in Eq. (|15p . 
We include only the range of coupling strengths where the 
fit is reasonable. It is noteworthy that our results with 
two lowest energies, Eq = — 4J and — 2J, show effective 
temperatures close to the degeneracy temperature, ap- 
proximately 2 J for this Hubbard system near half filling. 

We have not shown Toff for the highest energy we used, 
Eq = 1.77. This energy corresponds to the centre of the 
energy spectrum for all A plotted. At this energy, all 
the states of the subsystem have nearly equal statisti- 
cal weight at this energy. In other words, the effective 
temperature is nearly infinite. This is also reflected in 
the subsystem entropy (Fig. ^ which is close to log 16 
at Eq = 1.77, as expected for our 16-state subsystem at 
high temperatures. 

For systems exhibiting canonical thermalisation (small 
a^), the effective temperature Toff is approximately inde- 
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FIG. 8. Effective temperature T^s as a function of coupling 
strength A in the thermahsed regime for different composite 
energies Eq. Inset: example of fit of rss to the Boltzmann 
form 113. {U = J = 1, 5b = 0.5, L = 9.) 



pendent of the coupling strength up to A ~ 1. In fact, 
this effective temperature is close to the canonical tem- 
perature defined in Eq. (O by counting bath states in 
the limit of A — )■ 0, as reported in [s^. As already men- 
tioned, we find in this regime that the subsystem entropy 
(Fig. ini is also roughly independent of A. 

We will now turn to the crossover from non- 
thermalisation to thcrmalisation as we increase the cou- 
pling strength from zero. We can choose a rough measure 
of the threshold, Ath, for this crossover as the coupling 
strength at which a^^ drops below 25%. Alternatively, 
we can use the coupling strength at which the subsystem 
entropy reaches a plateau in Fig. [51 At £'o = —2 J, we 
find Ath — 0.05. At the lower energy Eq — — 4J, Ath is 
higher at approximately 0.1. At the energy Eq — 1.77 J 
corresponding to the centre of the spectrum, Ath is small- 
est at 0.03. The crossover between memory and lack of 
memory of the initial state also occurs around this char- 
acteristic coupling Ath- (We discuss this criterion fur- 
ther in Section |Vl) That thcrmalisation does not occur 
for small coupling strengths is because of the finite level 
spacing, A, in the finite-size bath. Physical intuition 
might suggest that subsystem-bath couplings, however 
weak, allow relaxation in subsystems. This is a reason- 
able assertion for systems with macroscopic baths where 
the bath spectrum is quasi-continuous. However, for a 
small system with a non-zero level spacing at weak cou- 
pling, the eigenstates of the composite system may only 
be slightly perturbed from the decoupled subsystem-bath 
product states \sb) if the typical matrix elements mix- 
ing these product states are small: {sb\XV\s'b') <^ A. 
In this weak-coupling limit, thcrmalisation cannot oc- 
cur from an initial subsystem state |0)sj when the com- 
posite eigenstates are all close to product states of the 
form 10)5 (8) |6)_B- The system would retain strong mem- 
ory of the initial state. Therefore, we expect a non-zero 



threshold for thcrmalisation for a finite system. We will 
examine more quantitatively the overlap of the compos- 
ite eigenstates with the decoupled product states in Sec- 
tion IIVI and we will compare the empirical Ath extracted 
here with a theoretical estimate in Section llVBl 

Let us now turn to the strong-coupling regime of 
A 1. As already discussed in Section UlI Al the sys- 
tem does not reach a steady state at very high A, and so 
it is not thermalised. We believe that this is a boundary 
effect in the sense that the dynamics in our 'subsystem' 
consisting of sites 1 and 2 become altered at very large 
A because of the very large hopping on the links between 
sites 2 and 3 and between sites 1 and L. As already 
discussed in Section Hi Al single-particle states localised 
on these links become visible as a feature the composite 
density of states at A = 10 (Fig. [3]). We believe that 
the four sites (i = L,l,2,3) will thermalise as a cluster in 
the sense that it has a canonical reduced density matrix, 
provided that the bath of size L — 4 is sufficiently large. 
Nevertheless, since the eigenstates of the two-site cluster 
and the four-site cluster are very different at large A, the 
thcrmalisation of the four-site cluster does not imply a 
diagonal RDM for the two-site cluster. In any case, we 
wish to make the point that this lack of thcrmalisation 
at large coupling is qualitatively different in origin from 
the lack of thcrmalisation at small coupling. 

It is interesting to examine more closely the depar- 
ture from thcrmalisation as we increase A in the range 
of A between 1 and 3 ioi U — J (see a^: in Fig. ^ . In 
this crossover region, we find steady states that have lost 
memory of the initial state (small Ar) but these states 
deviate from the canonical thermal state w, as can be 
seen in a rising a^^ as A is increased beyond unity. More- 
over, the RDM has a reasonable fit to the Boltzmann 
form (fT5|) . although the fitted temperature departs sig- 
nificantly from the canonical temperature ([7]). One can 
say that the system is still in an 'effective' thermal state 
in this crossover regime. We will return to this in Section 

ImcTI 

Interestingly, we observe that this crossover regime 
tracks closely a decrease in the density of states of the 
composite system. The density of states (Fig. |3]) can be 
approximated as a Gaussian: 



giEo) (X exp 



(£^0 — Eq)^ 



(16) 



where Eq is the energy of the band centre, and asw can 
be used as a measure of the width of the Gaussian. We 
see in Fig. [HI that (Tbw rises sharply as we increase A be- 
yond unity, similar to the behaviour of the fitted effective 
temperature (Fig. |6]). In fact, Teff cx crgy^,{X), as seen in 
Fig. ini where we show the two quantities normalised to 
their (A-independent) values at small A. Note that, at 
weak coupling and at fixed energy Eq, the derivative 



dloggiE) 



dE 



E — En 



En — En 



' BW 



(17) 
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FIG. 9. Gaussian width of the density of states, asw (sohd 
circles) and effective temperature Toff (hoUow squares) as a 
function of coupling strength A. {Eo — 0, U — J — 1, 5b ~ 
0.5, L — 9.) Both quantities are normalised to their values at 

Ath- 

can be associated with the inverse temperature of a bath 
of size L = 9. In other words, it appears that the effective 
temperature of the subsystem is better described by the 
canonical temperature of the whole system, instead of 
just the bath. This result is not surprising in this regime 
where the coupling of our 2-site subsystem to the L = 
7 chain is of order unity, since the distinction between 
subsystem and bath is blurred. 



D. Dependence on Interaction Strength 

We now turn to the effects of the particle-particle in- 
teraction strength, U, on thermalisation. In the results 
which follow, the coupling strength is fixed, as previ- 
ously, at A ~ 0.5 and we will also fix the bath width 
at Sb — 0.5J. The energy of the composite system, Eq, 
will be fixed such that it is always at the peak in the 
centre of the composite spectrum, at Eq « 2U. This 
is necessary since the shape of the spectrum is a strong 
function of U and the density of states at a given en- 
ergy can vary significantly. The effects of the interaction 
strength on the composite density of states are shown in 
Fig. m We find that for U > 4 J, peaks separated by U 
appear. If comparisons were to be made between differ- 
ent U for initial states at fixed Eg, the features in the 
density of states which evolve with U would introduce 
unwanted artefacts. Even in the centre of the spectrum, 
it should be noted that there is a fall in the density of 
states at the central peak at Eq « 2U, which occurs over 
a range 1 < U/ J < A due to an overall broadening of the 
density of states (see Fig. [3]). 

To measure thermalisation at different U, we will again 
employ the measures Ar and a^, as defined in Eqs. 
and (1121). The effective temperature is not shown since 
the initial-state energy is at a maximum in the density 
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FIG. 10. Memory of initial state, Ar , and closeness to the 
canonical thermal state, a^^, as a function of interaction U/J 
for A — 0.5. The composite energy Eq is chosen to be fixed 
on the central maximum in g{Eo), which lies close to 2U. 
Sb = 0.5J. 



of states which corresponds to infinite subsystem effec- 
tive temperature. In Fig. [TUl we demonstrate that ther- 
malisation, independent of the initial state, is found for 
U > 0.1. There is a broad minimum in plots of both Ar 
and , defined by the lack of thermalisation at small U 
and a small increase in the plotted quantities over the 
range 1 < C//J < 5. 

The slight increase in Ar and a^j above U c^i J coin- 
cides with the falling density of states in the centre of the 
spectrum shown in Fig.|3l So, the increase may be partly 
associated with the reduction in the number of states in 
the fixed bath window of our initial state. 

The behaviour at small U < O.IJ cannot be similarly 
related to the density of states. However when U — Q, 
the nature of the coupling is very different because the 
bath states are Slater determinants single-particle states. 
The single-particle level spacing is large compared to A J 
if A <C 1. So, we expect that thermalisation is poor for 
small non-interacting systems. In other words, for small 
[/, we need larger system sizes to observe thermalisation. 
We will present our data for different system sizes in the 
next subsection (Fig. [T^ . 



E. System Size Dependence 

It is interesting to study thermalisation as a function of 
system size. Owing to the exponential dependence of the 
Hilbert-space dimension on lattice size, it is not possible 
to find the full spectrum of large lattice. We will instead 
concentrate on the loss of thermalisation as we reduce the 
system size. If we use even smaller systems, reducing the 
number of sites rapidly leads to Hilbert spaces so small 
that thermalisation is not observed at all. Nevertheless, 
our results show that thermalisation is possible in sur- 
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prisingly small systems, as long as i7 ~ J so that the 
system is not close to the non-interacting limit, and as 
long as the density of states is not too low. To attempt 
to see some effects of reducing system size on thermalisa- 
tion, we consider composite states prepared with energies 
Eq in the centre of the band where the density of states is 
highest. As with our studies of the effects of interaction 
strength, this also eliminates unwanted effects due to the 
changing bandwidth with system size. 




FIG. 11. Closeness to canonical thermal state, aui, as a func- 
tion of coupling strength A for different system sizes L, for 
composite energies _Eo in the centre of the composite en- 
ergy spectrum. The number of particles was selected to keep 
= with the number of particles equal to L and L — 1 
respectively for even and odd L. {U = J = 1, 5b ~ 0.5.) 



First of all, let us consider how our results in Section 
nil CI for the dependence on coupling strength changes 
with system size. Shown in Fig. [TT] are plots of cr^ for 
different lattice sizes down to six sites. In each case, the 
subsystem size was fixed at two sites and the initial bath 
width was fixed at Sb — 0.5 J. Interestingly, for A ~ 1, 
thermalisation is maintained down to a four-site bath. 
However, the range of couplings over which thermalisa- 
tion occurs is greatly reduced. We return to system-size 
scaling in section Ivl 

We can also see how our results in Section IIII Dl for 
the dependence on interaction strength, U, change with 
system size. We see in Fig. [12] that the larger systems 
have a wider range of interaction strengths over which 
the system approaches the canonical thermal state (small 
aui)- Moreover, auj is lower for larger systems at a given 
U . This is consistent with our expectation that weakly- 
interacting systems require larger system sizes for ther- 
malisation. We will explore system-size scaling in Sec- 
tion El 




FIG. 12. Closeness to canonical thermal state, cT u: , as a 
function of interaction U/ J for different system sizes. A = 0.5. 



F. Dependence on Initial Bath State 



The thermal state should not depend on the micro- 
scopic details of the initial bath state. We will now 
demonstrate that the thermalisation behaviour found at 
long times is independent of the initial bath state. More 
specifically, we will vary the energy width (5_b of the initial 
bath state. For all of the numerical results presented thus 
far, we have considered initial states of the form (|4]) where 
the initial bath state is a pure state, with components in 
the bath eigenbasis non-zero only in a window of width 
Sb = 0.5 J. This was chosen because it is small compared 
with changes in the density of states. In Fig.fTSl we show 
plots of (T(j, Ar and Toff against A for values of (5b span- 
ning almost two orders of magnitude. In other words, 
these are results for vastly different bath states with the 
only constraint that they should be centered at the same 
energy. 

We find that the thermalisation behaviour is essentially 
(^B-independent for a broad range in Sb- Remarkably 
even up to 5b = 8 J, approximately half of the width of 
the composite eigenspectrum, we see Sb makes virtually 
no difference to the initial-state memory, quantified by 
Ar, and the effective temperature Toff. The distance to 
the thermal state at long times is modified slightly by 
choosing a very large Sb, but it should be noted that uj 
is itself dependent on the energy width of the state when 
this becomes large on the scale of changes in the density 
of states. 

Conversely, we can make Sb so small that there is just 
one initial bath eigenstate in the initial product state 
and virtually identical behaviour to Fig. [T2]is seen when 
A > 1. However, for smaller values of A, fluctuations 
appear as a function of A, thus necessitating a finite Sb- 



11 



1.2 D 

1 

0.8 


^ ' ■ ■ 5b = 6'.'i - 

\ 5o = 0.5 - 
" \ Sb = 2 - 

- \ Sb = 8 - 


— B 

A 

— V 

- 


< 0.6 






0.4 






0.2 



1 — 1 — 1 — v'Srai-ei 1 mwa b 





0.01 0.1 1 

X 



0.8 
0.7 \ 
0.6 
0.5 


, ' ■ ■ 5b = 6:'i - 

5b = 0.5 - 

\, Sb = 2 - 
- \ Sb = 8 - 


— B 

A 

^ 

— e — 


0.4 






0.3 






0.2 






0.1 









1 





0.01 0.1 1 

X 




X 



FIG. 13. Ar, and T^s as functions of coupling strength A 
for different bath- window widths 5b- The composite energy 
Eo = ~2J. The average level spacing A/J 10"^. (J = 1) 



G. Eigenstate Thermalisation 

We now discuss our results in relation to the eigen- 
state thermalisation hypothesis (ETH). As discussed in 
Section Hi CI this hypothesis requires the eigenstate ex- 
pectation values of the subsystem projection operator, 
Ps, (defined in Eq. (dU])) to depend only weakly on the 
exact choice of the eigenstate \A). 

We expect ETH to be valid in the regime where we 
found thermalisation in the previous sections — for U = 
J, this regime covers a wide range of coupling strengths, 
0.1 < A < 2, with A ~ 1 exhibiting behaviour closest to 



the canonical picture of thermalisation. So, we will now 
study the dependence of eigenstate projections 
on the coupling strength at U — J. We will focus on the 
projection on to the ground state of the subsystem in the 
{hs = 2, = 0) sector at ?7 = J. 

'Perfect eigenstate thermalisation' corresponds to the 
projection values forming a smooth quasi-continuous 
function of composite eigenenergy Ea- When this occurs, 
complete independence of the initial subsystem state ex- 
ists. Fig. UM shows histograms of (A|Ps=i|A). We see 
that there is some scatter in (A|Pi|A) for different eigen- 
states \A) that are close together in energy. There is the 
least scatter when the subsystem is closest to the canon- 
ical thermal state (small Ar and a^^) at A ~ 1 for U — J. 
Greater scatter in the values of (A|Pi|A) is found when 
the system starts to lose thermalisation (by our other 
measures of thermalisation), at small A < Ath = 0.1 and 
at large A > 2. 




-6 5 10 15 -5 5 10 15 

FIG. 14. Histograms (plotted as colour scale) of eigen- 
state projections pi — {A\Pi\A) on to the two-site subsys- 
tem ground state, for different subsystem-bath coupling A. 
{U = J = 1, Sb = 0.5J.) 

Let us examine the case oi X = 0.5 at U = J more 
closely. This shows little scatter, and hence good eigen- 
state thermalisation, over a wide range of energies. We 
quantify the extent to which eigenstate thermalisation 
holds by measuring the mean and the standard devia- 
tion, ctep , of the scattered values in each vertical column 
of histogram bins on the plot in Fig.[T31 This is computed 
using values within an energy window of 0.5 J. (This re- 
duces the fluctuations in the measured (Tep-) Fig. [TSl 
shows the positions of the mean values of (v4|Pi|^) and 
the positions of ±(Tep from the mean. We see that there 
is indeed good agreement between the mean eigenstate 
projections at any particular energy Eq and the canoni- 
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cal thermal value wn as defined in ([6]). 




FIG. 15. The same histogram of eigenstate projections pi — 
{A\Pi\A), from Fig. [14] with A = 0.5 is shown with a mod- 
ified colour scale, comparing the thermal values lju (black 
line) with the mean eigenstate projection from each vertical 
array of histogram bins (green line). Positions of one stan- 
dard deviation (ctep) either side of the mean of the eigenstate 
projection are also shown (blue lines). 

Let us now study the departure from eigenstate ther- 
malisation, measuring it by the increase in the scatter in 
the eigenstate projection values. To reduce any bias due 
to changes in the density of composite states with A, we 
used eigenstates at the energy Eg = 1.77 J which is near 
the maximum in the density of states for all coupling 
strengths considered here. The results are presented in 
Fig. 1161 First of all, we observe that a minimum in ctep 
indeed occurs over the same range of coupling strengths 
where other measures of thermalisation also show that 
the subsystem is close to a canonical thermal state. Sec- 
ondly, we find that, as the subsystem departs from eigen- 
state thermalisation at low coupling strengths, the in- 
crease in (Tep with decreasing A obeys the relationship 

CTEP oc i for A < 1. (18) 
A 

We will discuss this scaling in Section IIVBI 

We also lose eigenstate thermalisation if we increase 
the coupling strength to A ^ 1. As discussed before, we 
believe that this is a particular feature of our model where 
the properties of the coupling dominate the Hamiltonian. 

Finally, recall that we found in Section llll CI that. as A 
is increased beyond unity at U — J, the steady state 
of the subsystem departs from the canonical thermal 
state but the RDM follows a good fit to the Boltzmann 
form. This seems to indicate that the subsystem is in 
an effective thermal state that is non-canonical. We can 
see an indication of this crossover regime in Fig. [14] for 
1 < A < 3, where the scatter in the eigenstate projections 
is still relatively low, but the mean eigenstate projections 
as a function of Eq depart significantly from the canoni- 



cal thermal value ujn, in contrast to the case at A = 0.5 
(Fig. [13. 

In summary, we have shown that eigenstate thermal- 
isation holds and agrees well with other measures of 
thermalisation. We will demonstrate in Section IIVBI 
that the statistical behaviour of the eigenstate projec- 
tions (AlPsI^) is consistent with a simple model of the 
eigenstates \A) as random vectors in the basis of the 
subsystem-bath product states \sb). 




0.1 1 10 



FIG. 16. The spread of projection values, (tep, as a function 
of subsystem-bath coupling strength A (hollow squares). The 
values of (tep were found by averaging over composite energies 
in a window of width 0.5J centered on _Eo ~ 1.77. Solid line: 
theoretical estimate (|27[) for small A illustrating the scaling 
Ctep oc A~^. 



IV. EIGENSTATE OVERLAPS AND THE 
EIGENSTATE THERMALISATION HYPOTHESIS 

In Section IIII Gl we showed that the eigenstate ther- 
malisation hypothesis holds over a wide range of pa- 
rameters for our Hubbard-model system. Deutsch [l7[ 
and Srednicki [ll] have suggested that eigenstate ther- 
malisation occurs if the composite system is quantum 
chaotic. This was demonstrated theoretically for weak 
subsystem-bath coupling using results for the eigenstates 
of generic (random) Hamiltonians. In this section, we 
summarise these arguments and demonstrate numerically 
their agreement with our results for our Hubbard-model 
lattice. We will then use this framework to explain the 
dependence of eigenstate thermalisaton on the coupling 
strength discussed in section llll Gl namely the scaling of 
the spread of eigenstate projections, ctep, with coupling 
strength A for A < 1. 
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A. The Overlap Distribution and the Local 
Density of States 

To be more specific, eigenstate thermalisation is con- 
cerned with the eigenstate projections (A|Ps|A) = 
X]j,(A|s6)(s6|^). So, we need to understand the over- 
lap {sb\A) of the eigenstates |^) of the composite system 
at non-zero couphng with the eigenstates \sb) of the de- 
coupled system at A = 0. As for the eigenstate projec- 
tions the overlaps will fluctuate if we change 
Eji or Esb- However, we can study averages over energy 
windows that are narrow on the scale of variation in the 
density of states but contain enough states to smooth out 
fast fluctuations. 

The overlaps themselves are not invariant under a 
global gauge transformation and so should have mean 
zero. Let us consider first the squared overlap |(A|s6)p 
whose average is the variance of the overlaps. This can 
be interpreted as the weight of the product state \sb) at 
energy Est, in the decomposition of the eigenstate \A) at 
energy Ea using all the product states as the basis. In 
this picture of the composite eigenstate in energy space, 
|(^|s6)p is called the 'local density of states'. 

We now discuss some known results for the local den- 
sity of states. For a coupling XV between subsystem and 
bath, we expect that an eigenstate |^(A)) at Ea will con- 
sist mainly of product states \sb) with energies Esb close 
to Ea- If the coupling is not strong (A < 1), the energy 
range should scale with the strength of the coupling ma- 
trix elements A|(sy!iSB|y|s6)| where \sAbA) is the product 
state corresponding to |^) in the limit A — t- 0. To lead- 
ing order in A, this can be written as A|(^|y|s6)|. In this 
weak-coupling regime, we expect that the density of state 
of the bath spectrum i s nearly constant over this range. 
Then, the mean value, |(A|s6)p, should be a strong func- 
tion of the energy difference AEAsb = Ea — Esb, but has 
only a weak dependence on Ea and Esb separately. For 
a generic random coupling, we expect a Lorentzian form 
in the dependence on the energy difference: 



alsb^\{A{X)\sbW 



y\{A\V\sbW 



wl 



{Ea ~ EsbY 



(19) 



WL = ^y9{Esb)mv\sbw 



where g is the density of states of the composite system 
evaluated at the total energy Esb, taken to be approxi- 
mately constant over the energy width of this Lorentzian 
so that g{EA) — g{Esb) in this range of energy. It may 
be related to straightforward perturbation theoretic re- 
sults Q to second order in A. This result was origi- 
nally established [l^ over half a century ago for a spe- 
cific model of random coupling. This result was later 
shown to hold more generally. However, we note that 
this result does not take into account the specific case 
of coupling a bipartite system. Moreover, Eq. ([19]) is 
only strictly accurate in the general case for energy differ- 
ences where \Ea — Eb\> Wl- At smaller energy scales, 
non-perturbative mixing occurs and it is no longer pos- 



sible to associate eigenstates with a specific unperturbed 
state. However, the presence of non-perturbative mix- 
ing between states separated by less than Wl leads us to 
make the assumption that structures in the coupling ma- 
trix, such as elements which are identically zero because 
of the precise nature of the Hubbard-model coupling, are 
washed out by this mixing. We should also stress that 
this Lorentzian form is not expected to hold when A ~ 1. 
However, this is sufficient for us to use this form in the 
discussion of this section where we are concerned with 
the behaviour of the system at weak coupling. We note 
that, at stronger coupling, we find a Gaussian form for 
the local density of states [1^ . 
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FIG. 17. Local density of states, (Txst, as a function of 
/S.E = Ea — Esb for the Hubbard model with U = J = 1, 
at weak coupling A = 0.1. The averaging uses the overlaps of 
all eigenstates \A), with |s)g = |t)t) and |&)s selected within 
energy J from the centre of the bath spectrum. Solid line: 
Lorentzian with width Wl as given by (|20[) . 

We can make an estimate for the width Wl in the 
Lorentzian form ([TO]) . As discussed above, the overlap 
(A|]/|s6) can be approximated by {sAbA\V\sb) to lead- 
ing order in A. So, we see that its mean square value 
should be the mean square value V"^ of an element of 
the coupling matrix {s'b'\V\sb) - We will see in Eq. ((3T|) 
in Section |V] that V"^ ~ JA/2 for an interacting system 
near half filling. If we further approximate the density 
of states g[Esb) with the average level spacing A, we see 
that 



Wl 



ttAV 



(for [/ - J) . 



(20) 



The local density of states, cr^s^' Hubbard model 

at A = 0.1 and C/ = J is shown in Fig. [T71 It fits well to 
the Lorentzian form (|19l) with Wl given by ([20)) . 

We can also discuss the full distribution of the overlaps 
{sb\A). For a system with time reversal symmetry, {sb\A) 
can be constructed to be real. Our numerical results for a 
system dX U = J and A = 0.1 are shown in Fig. [TH] This 
is a histogram using the overlaps of all the eigenstates \A) 
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FIG. 18. Two-dimensional histogram showing the distribu- 
tion of the overlaps {A\sb) as a function of the energy differ- 
ence Ea — Esb for the Hubbard model with U = J = 1, at 
weak coupling A = 0.1. The histogram includes the overlaps 
of all eigenstates \A), with \s)s = | tit) and |6)_g selected 
within energy J from the centre of the bath spectrum. The 
histogram bin widths are 0.02 and 0.002 on the energy and 
overlap axes respectively. 



with a subset of product states \sb) where the bath states 
are within an energy J of the centre of the bath spectrum. 
We can see that the width of the distribution is a strong 
function of the energy difference AEAsb between \A) and 
\sb). The widest distribution is found at Ea = Esb- In 
this case, the states \sb) effectively form a random basis 
for the eigenstates \A). 

The distributions appear to be controlled by a single 
variable, the local density of states. In other words. 



P{X = {sb\A)) = a^l,F{X/aAsb) 



(21) 



for a normalised distribution F{u) with unit variance. 
This is demonstrated in Fig. [19] for our data at five dif- 
ferent AEAsb- The data have been scaled using the ex- 
pected width aAsb given by (fT9|) and ([20| . So, this data 
collapse contains no adjustable parameters. The distri- 
bution F{u) has an excess kurtosis 7 = (w"*) — 3(?i^)^ ~ 2 
numerically. (7 would be zero for a Normal distribution.) 
In Fig. [HI we see that our data are well approximated 
by a hyperbolic secant distribution which has an excess 
kurtosis of 2: 



F{u) 



1 



2 cosh(7rM/2) 



(22) 



We should point out that the data collapse to this distri- 
bution fails at strong coupling. This may be due to the 
fact that the width of the distribution aAsb becomes large 
enough that each eigenstate \A) involves bath states in 
a wide range of energies over which the bath density of 
states varies significantly. 

We note that our form for the overlap distribution dif- 
fers from what may be expected from random matrix 



theory [l7|, [l8|, [23j for similar types of overlaps which 
suggests that they should follow a Normal distribution 
at weak coupling. This indicates that the overlap distri- 
bution may depend details of the coupling Hamiltonian 
or details of random matrix ensemble. 
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FIG. 19. Distribution of overlaps X — {A\sb} at weak cou- 
pling A = 0.1, scaled by the width aAsbi^E = Ea — Est) at 
different values of AE. (J — 1; same system parameters as in 
Fig. 1181 ) Solid line: hyperbolic secant distribution with zero 
mean and unit variance. Inset: Log plot of the same data. 



Accepting the distribution ([21]) as the distribution for 
the overlaps, the distribution P^^^ for the squared over- 
laps |(s5|v4)p can be derived: 



P(2) {\{A\sb)\^ ^Y) 



F{VY/aAsb) 

(JAsbVY 



(23) 



which has a mean of (JAsb ^^'^ ^ variance of (2 -I- 'y)o'Asb- 
In the case where the overlap distribution is so wide that 
\A) is effectively a random vector in the basis of \sb), we 
have a Porter-Thomas distribution for the local density 
of states. 

To summarise, we have shown that our numerics agree 
with results for the local density of states arising from 
generic random Hamiltonians. This controls the overlap 
distribution. We point out that in this simple picture 
of the statistics of the overlaps, any correlations between 
different eigenstate overlaps are implicitly neglected. We 
will now proceed to understand eigenstate thermalisa- 
tion in terms of this simple picture of eigenstate over- 
laps {ab\A)- 



Scaling of Eigenstate Thermalisation with 
Coupling Strength 



In Section ITlIGi we found that the degree of eigenstate 
thermalisation improves upon increasing the strength of 
the subsystem-bath coupling. In the previous section, we 
have seen (|19p that, in parallel, increasing the coupling 
strength broadens the local density of states. We expect 
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that, as the distribution of overlaps {sb\A) broadens such 
that more basis states \sb) participate in each eigenstate, 
the fluctuations in the projection values (A|Ps|A) will be 
reduced in accordance with the law of large numbers. We 
will show this to be the case and find the A-dependence 
for the spread of eigenstate projections ctep, found nu- 
merically in Section [ml 

From its definition ((TU| . the projection operator P, 
sums over all bath states. In our model of the overlaps, 
= |(s&|>l)p is a sum over many independent 
variables. Its mean and variance are given by 



2 



-(2 + 7)E 



4 

^Asb ■ 



(24) 



For the full distribution of these quantities, see Ap- 
pendix which applies the central limit theorem to our 
model distribution for the overlaps. 

First, we consider the mean hep, As for a given sub- 
system state s. From our model this is the sum 
over all bath states (with the necessary spin and parti- 
cle number) using a Lorentzian window of energy cen- 
tered at Cb = Ea — £s- In other words, the answer 
should be proportional to the bath density of states at 
Ea — £s- Furthermore, we have the normalisation condi- 
tion {A\A) = Y.sb ^\sb = J2s MEP,As = 1- So, at fixed A, 
/iEp should give the normalised probability of finding the 
subsystem in state s according to the Gibbs distribution: 



AiEP - E o-Asb - 
b 

gsiEA - es,N - ns,S 



b"Asb 



J2s J2b '^Asb 



sf) 



giEA,N, S^) 



uJss{Ea) (25) 



where Ea , N and are the energy, number and spin of 
the state \A), 5s(eb, rif,, s^) is the density of bath states 
with energy in an interval about e^, particles and spin 
sg. Thus, we see that a simple model of the eigenstate 
overlaps gives the canonical thermal distribution 

Next we address the spread of the projection values 
(TEp. Using the Lorentzian form (IT51) for the local den- 
sity of states with width Wl, and following the same 
approximations as above, the sum over a\ is: 



gBieb)Wl/7r^g^{EA) 
[Wl + {Ea - - eb)^Y 



deh 



gB{EA-£s) f {Wl/7r^)dei, 



gHEA) J [Wl + {EA-es-ebfV 

UJss{Ea) 



2^WLg{EA) 



(26) 



This means that the spread of the projection values is 
given by: 



\ {2 + -f)ujss{EA) _ 1 / (2 + 7)a.ep,a. 
""P-V 2.W,g{EA) --X^ Jg{EA) ^'^^ 



where we have used our estimate pOj) for the Lorentzian 
width Wl. Note that l/g{EA) is of the order of the 
average level spacing A. 

Therefore, we find that ctep is proportional to as 
shown numerically in Fig. 1161 Moreover, we see that 
(Tep ^ •\/A so that the fluctuations in the projection 
values are small for large systems, in accordance with 
the law of large numbers for a quantity that is a sum 
over many states. 

We stress that the above results hold for any dis- 
tribution of matrix elements (A|s6), of sensible form, 
where the central limit theorem applies. Furthermore, 
the result p5|) holds quite generally for any sensible 
form of (j\gb "which is a function of Ea — Esb with a 
peak at Ea — Egb = 0. Although we have not derived 
it here explicitly, it should also be noted that the off- 
diagonal elements of the reduced density matrix, which 
were found to be virtually zero numerically, are expected 
to be zero from this model of eigenstate overlaps. Indeed 
the mean values of the eigenstate expectation-values for 
off-diagonal elements are clearly zero due to the random 
sign of eigenstate overlaps. 



V. CROSSOVER TO THERMALISATION 

In this section, we will try to understand the on- 
set of thermalisation using simple theoretical arguments. 
In particular, we explore the effects of system size on 
the thermalisation for the cases of interacting and non- 
interacting fermions in the Hubbard model. The ex- 
tent to which the effects of system size on thermalisation 
may be seen numerically is limited, as discussed in Sec- 
tion nil El We proceed to identify a minimum coupling 
strength Anp below which thermalisation cannot occur. 

To begin to see thermalisation requires the coupling 
Hamiltonian XV to be big enough to mix the A = 
eigenstates non-perturbatively. We will then compare 
this theoretical estimate with our numerical results. 

For small coupling strength (A ^ 1), the overlap be- 
tween an eigenstate |^) and a subsystem-bath product 
state \sh) takes the form 



{sb\A) ~ SsAsSbAb + X 



+X' 



E 



{sb\V\sAbA) ' 
EsAbA ~ ^sb_ 

sb\V\s'b'){s'b'\V\sAbA) 



s'b'^SAbA 



SAbA 



Esb){E. 



SAbA 



Es'b') 



1 \{sb\V\sAbA)? 



2{Es 



EsbY 



(28) 



to second order in A, where the state |sa&a) is the com- 
posite eigenstate \A) to zeroth order in A. (Note that 
V has no diagonal elements in this basis.) The threshold 
for non-perturbative mixing may be considered to be met 
when the second-order term equals the first order term 
in magnitude. Note that the bath states |6)b coupled 
by V have different quantum numbers from the given 
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state \bA)B so that there is no level repulsion between 
\b)B and |6^)_b- We expect the nearest bath state is on 
average Ab/4 away in energy. Generically, this occurs 
around a coupling strength Anp which we define by 



Ann 



y2 



4 



(29) 



where As is the bath level spacing and V"^ = 
|(s'6'|y|s6)p is the typical magnitude of the square of 
a coupling matrix element. We will estimate these below 
for interacting and non-interacting systems. 

The quantity A,ip should be the coupling strength at 
which one starts to see a departure from complete mem- 
ory of the initial state at long times. We therefore expect 
that this quantity should be similar to the quantity Ath, 
introduced in Section IIIH which measures the crossover 
from the non-thermalised regime to thermalisation. Note 
that Ath has been defined with an arbitrary choice of a 
threshold for cr^ at 25%. Its actual value will change 
with the specific criterion chosen to mark this threshold. 
However, one can use the data from Fig. [TT]to show that 
the relative values for Ath for different system parame- 
ters are approximately the same for a range of choice of 
thresholds. So, it is reasonable to discuss a relationship 
between Anp and Ath- In particular, it is expected that 
the two quantities should be proportional to each other 
for a given subsystem size. 

The rest of this section is dedicated to understanding 
the scaling of Anp with system size. First we will con- 
sider the case of finite interactions U ^ J before, in the 
subsection following, discussing the case of virtually free 
fermions where U <^ J . 
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FIG. 20. Coupling matrix elements linking the = 3, = 
i and ns = 2, sf = 1 subsectors for the Hubbard model 
with U = J = 1. Colour scale indicates the magnitude of the 
matrix elements. The banded diagonal structure is typical for 
the coupling of all subsystem states. Sizeable matrix elements 
lie within a band of width 4J = 4, with some very small 
matrix elements lying outside of the band due to the finite 
interaction strength U . 



First of all, we estimate the magnitude of a matrix 
element {sh\V\s'b'). To keep the description straightfor- 
ward, we consider the case of exactly half filling. We 
will compute this from the average for the sum of all the 
squared matrix elements TiV"^ — J^ss'bb' \{sb\V\s'b')\'^ . 
The calculation of this trace can be found in Appendix 
m We find: 



A. System-Size Scaling for Interacting Fermions 



Tr(y2) ^ 2MJ2 



(30) 



We will now deduce the scaling of Anp with system size 
for the Hubbard model with interactions U ^ J. To find 
the theoretical scaling of Anp with system size requires 
a knowledge of the scaling of both the energy spacing 
between coupled states and the scaling of the magnitude 
of the typical matrix elements {sb\V\s'b') with system 
size. A characteristic submatrix of the coupling matrix 
{sb\V\s'b'), with s and s' fixed, is shown in Fig. [20l for 
the Hubbard model with interaction strength U = J. 
The non-zero elements of the coupling matrix form a 
band. This can be explained by the single-particle na- 
ture of the coupling. In the limit of zero interactions, the 
coupling involves a single particle hopping into, or from, 
one of the single-particle states in the bath. Therefore, 
the full width, 2W, for bath states into which a particle 
may hop is 4J, the single-particle bandwidth. The pres- 
ence of interactions preserves the banded structure of the 
coupling matrix and, provided U < J, the banded ma- 
trix is not significantly broadened beyond 4 J. However, 
when U ^ J, the details of single-particle bath states 
are blurred as the single-particle quasiparticle weight is 
significantly reduced from unity. 



where M is the dimension of the Hilbert space of the com- 
posite system. We now need to count the number of non- 
zero matrix elements in the coupling matrix {sb\V\s'b') . 
Since the coupling involves the hopping of a single par- 
ticle of a given spin state between the subsystem and 
bath, any given subsystem state s, will only have non- 
zero matrix elements with at most four other subsystem 
states s', corresponding to changing the particle number 
or spin by ±1. So, there should be approximately 4M5 
such non-zero blocks in the coupling matrix where Ms 
is the dimension of the subsystem Hilbert space. Each 
block has a banded structure similar to the one shown in 
Fig. [201 Note that the bath states, b and b' , connected 
by {sb\V\s'b') belong to subsectors of the bath spectrum 
with different quantum numbers. For a band of full width 
2W, the banded block should have Ms (211^/ As) non- 
zero elements where A^ is the average bath level spacing 
and Mb is the total number of bath states in the bath 
subsector of a given number and spin. So, the total num- 
ber of non-zero elements in the coupling matrix V is ap- 
proximately 4MsMB(2I^/Ai3) ~ 8MbW/A. Therefore, 
the mean squared value of each coupling matrix element. 
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is 



1/2 



2MJ2 



(31) 



using 2W = 4 J which is valid for the case in Fig. [20l 
where [/ ^ J as discussed above. So, from (|29p . we see 
that non-perturbative mixing occurs when A reaches 



A 



np 



Ab 
2J 



MgA 
2J 



(32) 



For the purposes of understanding how this threshold 
scales with system size, we have approximated the bath 
level spacing as a simple multiple of the average level 
spacing A of the composite system: A^ ~ Ms A. 

We have arrived at this condition using simple argu- 
ments based on perturbation theory. A similar criterion 
can be obtained using our results for eigenstate thermal- 
isation. In terms of the eigenstate projection values, the 
system does not thermalise if the spread of the projec- 
tion values, (Tep, becomes comparable to the mean /iEP- 
We have shown in (I^S)) that the latter gives the canon- 
ical state Ld which is of the order of 1/Ms where Ms is 
the number of states in the subsystem. Thus, from (P?]). 
the condition that ctep < Mep for a given total energy 
Eq and a subsystem state s becomes the criterion that 
A > XETniEo, Ss) where 



■\ETH 



2 + 7 



JgB{Eo 



(33) 



of the initial state which allows the onset of thermalisa- 
tion occurs at 



XnpiU - J) = 4V2i2- 



(36) 



Reassuringly, Anp tends to zero as L tends to infinity so 
that, for baths in the thermodynamic limit, arbitr arily 
small couplings lead to thermalisation, as we expect [41j . 
For a lattice of nine sites we estimate this sets Anp — 
0.054. 
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FIG. 21. A comparison of Ath, Anp and Aeth for different 
system sizes L when U — J. Ath is set as the coupling strength 
for which a^j falls to 0.25 for composite energies Eo in the 
centre of the band. Anp is found using ([32]) with (|34[) . Aeth 
is obtained from (1331). 



Since l/A^ is simply the average of gs over the bath 
sector, we see that the thresholds Anp and Aeth, based on 
different criteria, describe essentially the same crossover. 
^ETH is larger than Anp as might be expected since the 
latter marks the loss of memory of the initial state while 
the latter marks the onset of the canonical thermal state. 

It remains to establish how the level spacings A de- 
pends on the system size. Assuming the cosine dispersion 
for a tight-binding band and neglecting the broadening 
due to finite U, the many-body bandwidth for L sites 
is 8JL/tt at half filling. This is found approximately by 
finding the maximum and minimum composite energy 
eigenvalues, zLAJL/tt, by summing the energies of the 
L/2 highest and L/2 lowest single-particle eigenstates. 
Therefore, the mean level spacing is 



A = 



8JL 



Mtt 

At half filling, the Hilbert-space dimension is 

H 2 

M = 



(34) 



(35) 



Therefore, using p2p and Stirling's approximation for 
large L, we find that the threshold for the loss of memory 



We now compare Anp, for different lattice sizes, with 
Ath- To allow easy comparison between different system 
sizes, the composite energy in the centre of the band will 
be considered in each case. As already discussed, the 
value of aui where the value Ath is recorded is somewhat 
arbitrary. However, we find good agreement between Ath 
and Anp , as is shown in Fig. I21[ with Ath defined using a 
threshold of a^j = 25%. If thresholds other than cr^=25% 
are considered, we find that Ath changes approximately 
by a multiplicative constant for all L. Therefore, the 
good agreement between the explicit values for Ath find 
Anp is not a remarkable feature of Fig. [5T] However, 
that the two quantities are found to scale in virtually 
the same way for the limited numerical data available 
provides numerical evidence to support the system-size 
scaling of thermalisation (I55|) derived above. 



B. System-Size Dependence for Non-Interacting 
Fermions 



For the case of almost free fermions {U <C J), thermal- 
isation was not seen in the nine-site Hubbard ring. We 
now repeat the argument above for the case of negligible 
U. The major difference from the case of ?7 ~ J above 
is the structure of the coupling matrix. Shown in Fig. [51] 
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is the striped form of the couphng matrix. Without in- 
teractions, the bath eigenstates are simply Slater deter- 
minants of free-fermion single-particle states. Therefore, 
for each spin, the coupling Hamiltonian has non-zero ma- 
trix elements only at energies corresponding to the L — 2 
single-particle bath states for each spin. 
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FIG. 22. Coupling matrix elements linking the Ua — 3, Sg = ^ 
and ris = 2, = 1 subsectors for the Hubbard model with 
U — O.OIJ = 0.01. Colour scale indicates the magnitude of 
the matrix elements. The banded diagonal structure is typical 
for the coupling of all subsystem states. The finite matrix 
elements lie within a band of width 4J but, in contrast to 
Fig. 1201 the matrix appears striped. 

When considering the threshold for non-perturbative 
mixing, it must now be noted that coupled states differ 
not by the level spacing A, but by the bath single-particle 
level spacing Ai, where Ai « 4J/(i — 2). The magnitude 
of TrF^, as given by pop . is independent of U. Therefore, 
using ((32|) with A replaced by Ai, we see that Anp for 
U J should be given by 



A„p(C/« J) 



L-2 



(37) 



This yields a value Anp — 1.1 for the nine-site lattice. It is 
therefore clear why initial-state independence is not seen 
for the nine-site lattice when U J: the threshold for 
non-perturbative coupling occurs at a coupling strength 
some 20 times bigger than for the Hubbard model with 
interactions U ^ J. Eq. (P7)) indicates that the U J 
case needs a lattice with 2700 sites in order that Anp falls 
to the same value as for the C/ ~ J interacting case. 

We have argued that thermalisation for small systems 
occurs at smaller system sizes for the interacting system 
compared to the non-interacting system. We will now ask 
how the crossover from the U <^ J regime to the U ^ J 
regime occurs. This should occur when the width of the 
stripes in the coupling matrix elements seen in Fig. [25] 
becomes comparable to the single-particle level spacing 
Ai. The stripe width should be the quasiparticle decay 



rate due to interparticle collisions. At small C/, the decay 
rate can be estimated using Fermi's golden rule. The 
matrix elements are proportional to U and the density 
of single-particle states is proportional to 1/J. So, for 
single-particle energies far from the Fermi level so that 
we can neglect effects from Pauli exclusion, the decay 
rate should be ~ U"^ / J . This becomes comparable to Ai, 
when U reaches C/th ^ VTAi cx 1/ \JL — 2. This estimate 
give the scale for the interaction strength beyond which 
we see thermalisation in the numerical results shown in 
Fig. [T^ Wc have only four system sizes and there is 
a strong even-odd effect in the system size, making it 
difficult to verify our prediction quantitatively. 



VI. EXPERIMENTAL IMPLICATIONS 

Models such as the Hubbard model studied in this work 
can be simulated readily using cold atoms in optical lat- 
tices [13 . Thanks to recent rapid progre ss in addressing 
single sites in optical lattices [20|, I43l - l45| , models similar 
to those we studied here can now in principle be im- 
plemented and measured in systems of ultracold atoms 
trapped in optical lattices. In particular, single-site imag- 
ing capability means that atom occupation (albeit up to 
number modulo two) and the spin species can be de- 
termined accurately at a few lattice sites that will form 
the subsystem. This means the state of the subsystem 
can be probed directly. Single-site addressability means 
the subsystem and the bath can be initialised with pure 
quantum states with well-defined number and spin. Fur- 
thermore, instead of focusing the probe laser beam on a 
single site, the laser can be aimed accurately (to within a 
tenth of the lattice spacing Q) between two neighbour- 
ing sites, to tune the lattice potential locally and thus 
adjust the coupling A J between the subsystem and bath. 
Both A > 1 or A < 1 regimes can be accessed with a blue- 
or red-detuned laser focussed between the sites. 

We expect our findings to be seen for bath state with 
a relatively well-defined energy (which overlaps bath 
eigenstates only within a range of energies much smaller 
than the many-body bandwidth), far from a strongly- 
correlated ground state. We should point out that, al- 
though we focussed on a specific Hamiltonian with spe- 
cific initial conditions, we believe that our results are gen- 
erally applicable. For instance, we obtain similar results 
for Bose and fermion Hubbard models. Also, although 
we have used an initial bath state (U) consisting of bath 
eigenstates within a narrow energy window, our results 
(Fig. [T3l) are not sensitive to the width of this energy 
window. Our results should hold for initial states span- 
ning a larger window of bath energies which would be 
easier to prepare experimentally. Our results also do not 
change if we introduce random coefficients in the linear 
superposition of bath eigenstates or change the shape of 
the window, consistent with the results of Ref jld^. 

Moreover, experimental systems will be larger and con- 
tain many more atoms than in our simulation. We have 
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shown that the threshold Ath in the couphng strength is 
exponentiahy smah in the system size for large systems. 
So, we believe it would be possible to see thermalisation 
at smaller A in experimental systems, even if the initial 
bath state is simply a single eigenstate. 

We should also ensure that the time needed for ther- 
malisation to be seen should be within the lifetime of 
an optical lattice experiment, typically hundreds of mil- 
liseconds or more. Our previous work [36j for systems 
with [/ ~ J shows that the relaxation towards equilib- 
rium occurs with a relaxation rate ^ for weak cou- 
pling (showing exponential decay) and ~ A J when A ^ 1 
(showing Gaussian decay). (See Fig. 4 of [l!].) For cur- 
rent optical lattice experiments with '^''K, using an optical 
lattice laser wavelength of 1064 nm, the hopping matrix 
element J is approximately 380 Hz for a laser strength 
Vb of 5 times the recoil energy {En), and ~ 100 Hz for 
Vb — lO-Efl. Hence, expected relaxation time scales when 
Vb = bEfi range from about 3 ms at A = 1 in the Gaus- 
sian regime, to about 30 ms at A = 0.1 in the exponential 
regime. (The corresponding time scales for Vb = lOi^i? 
range from 10 ms to 100 ms.) For the other commonly 
used species ^Li, the lighter mass means shorter time 
scales than for '^^K. For an optical lattice laser wavelength 
of 1064 nm, the corresponding Gaussian relaxation time 
scale is about 0.4 ms and 1.5 ms for exponential regime, 
with Vb — 5Ei^. Hence, for optical lattice laser strengths 
that are not too large, the relaxation times are well within 
experimental lifetime of the cold atom systems. 

Finally, we point out that it is not necessary to use a 
ring or one-dimensional geometry (as studied in this pa- 
per) to see the thermalisation physics we have presented. 
As long as the inelastic scattering length is small com- 
pared to the size of the bath, we expect the qualitative 
aspects of thermalisation to survive, although the pre- 
cise values for the various crossovers and thresholds will 
change depending on the dimensionality of the bath and 
the nature of subsystem-bath coupling. 



VII. CONCLUSIONS 

We have presented an account of the thermalisation of 
a local subsystem within a closed quantum system de- 
scribed by a lattice of interacting fermions. The subsys- 
tem thermalises in the sense that its reduced density ma- 
trix approaches the form expected for a canonical ther- 
mal ensemble. This thermalisation occurs over a wide 
range of system parameters for surprisingly small sys- 
tems. The equilibrium state depends very little on the 
strength of subsystem-bath coupling provided it meets 
the following two conditions. Most importantly, the cou- 
pling strength needs to be large enough to mix the eigen- 
states of the uncoupled system non-perturbatively. Sec- 
ondly, the coupling strength must not be so large that the 
boundary effects associated with the coupling dominate 
the behaviour. We also find that small lattice clusters 
thermalise for a range of interaction strengths, provided 



U is large enough that the system is away from integra- 
bility at [/ = 0. We were also able to demonstrate that 
the energy width of the initial pure state of the bath has 
virtually no effect on the subsystem state at long times. 
This was found for a range of energy widths spanning 
nearly two orders of magnitude. 

Numerically, we demonstrated the relationship be- 
tween subsystem thermalisation at long times and the 
eigenstate thermalisation hypothesis. We further quanti- 
fied the extent to which eigenstate thermalisation holds 
by measuring the spread of eigenstate expectation values 
for subsystem occupation probabilities. Using generic re- 
sults for the eigenvectors of perturbed quantum systems 
in random matrix theory, we were able to derive theoret- 
ically a coupling-strength threshold Aeth for thermalisa- 
tion which is in qualitative agreement with our numerical 
threshold Ath- This establishes a link between the eigen- 
states of weakly-coupled bipartite quantum systems and 
eigenstate thermalisation. As this result employed only 
random matrix theory, our conclusions should be quite 
general for non-integrable systems, provided that the sys- 
tem is prepared at an energy far from the ground state 
where correlations may become important. 

We were also able to understand the system-size scal- 
ing of the breakdown of thermalisation seen in our nu- 
merics for interacting fermions by considering a coupling- 
strength threshold, Anp, below which non-perturbative 
mixing of A = eigenstates does not occur. We demon- 
strated that this non-perturbative threshold Anp has vir- 
tually the form as Aeth- Moreover, these have the same 
system-size scaling as the empirical Ath- 

We deduced that these thresholds for thermalisation 
should tend to zero exponentially in the system size. 
We also attribute the lack of non-perturbative mixing 
as the reason for the lack of thermalisation for the 
weak-interaction limit of the small systems we stud- 
ied. For very large Hubbard rings, we predict that non- 
perturbative mixing does occur for any non-zero interac- 
tion U. 

During preparation of this manuscript we became 
aware of unpublished work by Neuenhahn and Mar- 
quardt 46] which also studies eigenstate thermalisation 
using random-matrix-theory results for eigenstate over- 
laps. The authors study the momentum distribution of 
interacting fermions on an entire closed system, in con- 
trast to the local observables on bipartite quantum sys- 
tems considered in this work. 
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Appendix A: Overlaps and related distributions 

We start with the distribution of overlaps X ~ {sh\A) 
in Eq. ((2T|) with zero mean, variance {X"^) — (^\sb ^^"^ 
fourth moment (^'*) = (3 + ^)<j\g-)j- Since the projec- 
tion is given by ^j^{A\sb) {sh\A) , we first consider the 
distribution, P'^) foj. \ l^A\sb)\'^ : 

P^^'>{\{A\sh)\^ =YAsb) 

/OC 
dXAsbP{XAsb)5{YAsb - ^LJ 
-oo 

.,F\^^^\ (Af) 



f 



(^Asb^/YAsb \ ^Asb 

which has mean cr^^j and variance {2-\-^)a\^^^. Then, the 
eigenstate projection values have the distribution P^^, 
given by 



M. 



J " B / poo \ j'-f-B 

n / dY,P^'\Y,)) 5{w-J2y,) (A2) 

1=1 ^-^O ^ fc=l 

where Mb is the total number of bath states and, since 
A and s are fixed, the notation is abbreviated such that 
YAsb — ^ ^fc- Since this is a sum of many independent 
random variables, albeit from different probability distri- 
butions, it is reasonable to ask if a central limit exists. 
Indeed, the Lyapunov condition for a generalised cen- 
tral limit does hold [131 • To find this central limit, we 
adopt the standard procedure of factorising the integrals 
in Fourier space. Upon taking the Fourier transform 



dWe^'^P^^iW) 



(A3) 



the contributions from each of the bath states factorise 
such that 



(A4) 



where 



Pf (fc) 



1 + i/ccrisb - 



(A5) 



where the series has been truncated to second order in k. 
The logarithm of P^^(fc) takes the form of the series 



logP^P(fc) 



b 



(A6) 



The coefficient to the term linear in k is simply i/^ep and 
the coefficient to the k^ term is — (Tgp/2 where /iep and 
(Tep are defined in (|24l) . We have dropped terms of higher 



order of the form X]6(^'''Asfc)"- Using and following 
the same argument that leads to (pS)) . 



fc"MEP,As 

[Wl9{Ea)]^- 



(A7) 



where g{EA) is the density of states at Ea and 0Jss{Ea) 
is the reduced density matrix for the canonical thermal 
state ([6]). Therefore, we see that the truncation of the 
series is reasonable for k <^ gWh oc A^J/A. 

Upon re-exponentiating the series, we see the bulk 
of the distribution P^^{k) may be described accurately 
with A: up to the scale of I/ctep oc ( J/ A) since the cen- 
tral limit only breaks down at fc ~ J/ A. (This condition 
is readily met in our numerics when the coupling strength 
A is large enough for the subsystem to approach thermali- 
sation.) Exponentiating and inverting the Fourier trans- 
form yields the distribution for eigenstate expectation- 
values: 



1 



V27R7eP 



exp 



{W - MEP) 

2a|p 



, (A8) 



which is a Normal distribution with mean /iep and vari- 
ance ctep- 



Appendix B: Coupling matrix 

In this section, we estimate the magnitude of a matrix 
element of the coupling matrix V as defined in ([3]). As 
discussed in Section |Vl the coupling matrix involves only 
single-particle hopping between the subsystem and the 
bath. So, it should connect states not further apart in 
energy than the single-particle bandwidth 4 J. 

We will consider the coupling matrix to be a banded 
matrix where the non-zero elements form a band of full 
width 2W = 4J. While enumerating the size of individ- 
ual matrix elements is not possible without full diagonal- 
isation of the A = Hamiltonian, the quantity Tr V'^ is 
basis-independent and may be found readily in the Fock 
basis, with states |Pi), where particles are localised. In 
this case: T^{V^) =Y..A{m\F,)? ■ 



A 

Subsystem 
site i=2 Q 

D 



4 



Bath 
site i=3 



4 f 



FIG. 23. A diagram showing four possible occupations of 
two sites (across the coupling link at i=2 and 3) by spin- 
up fermions, irrespective of the configuration of spin-down 
fermions on these sites. 
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The matrix V does not change the total particle num- 
ber. To keep the description straightforward, we con- 
sider the case of exactly half filling. This is demon- 
strated in Fig. [231 For each basis state {Fi \ , there are 
at most only four other basis states \Fj) which are re- 
lated by hopping a single fermion (spin up or down) be- 
tween the subsystem and the bath via either one of the 
two subsystem-bath links. As the lattice is taken to be 
exactly half-filled, for each spin and for each topological 
link between subsystem and bath, half of the Fock states 
have a filled site adjacent to an empty site across each 
coupling link. This diagram shows four possible occupa- 
tions of two sites (across the coupling link at i—2 and 
3) by spin-up fermions, irrespective of the configuration 



of spin-down fermions on these sites. At half filling, the 
full L-site Fock states may be divided up into four groups 
containing equal numbers of states, with each group hav- 
ing the spin-up occupations A, B, C and D (as labelled 
in the figure). Each state in groups A and B can cou- 
ple to one other Fock state, with matrix element A J, but 
states in groups C and D couple to no other Fock states. 
Spin-down fermions do not affect these matrix elements. 

Therefore, each spin and each subsystem-bath link con- 
tributes MJ^/2 to the trace where M is the dimension 
of the Hilbert space of the composite system. There are 
contributions from two links and two spin species. Hence, 
we obtain 

Tt{V^) = 2MJ^ . (Bl) 
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